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Abstract 

We consider the clustering of Lennard- Jones particles by using an energetic connectivity criterion 
proposed long ago by T.L. Hill [J. Chem. Phys. 32, 617 (1955)] for the bond between pairs 
of particles. The criterion establishes that two particles are bonded (directly connected) if their 
relative kinetic energy is less than minus their relative potential energy. Thus, in general, it depends 
on the direction as well as on the magnitude of the velocities and positions of the particles. An 
integral equation for the pair connectedness function, proposed by two of the authors [Phys Rev. E 
61, R6067 (2000)], is solved for this criterion and the results are compared with those obtained from 
molecular dynamics simulations and from a connectedness Percus-Yevick like integral equation for 
a velocity-averaged version of Hill's energetic criterion. 
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I. INTRODUCTION 



The concepts of clustering and percolation have been widely used in order to explain 
several phenomena in very diverse areas including Physics, Chemistry, Biology, Geology, So- 
ciology and Economics. In particular, with reference to chemical-physics, phenomena such as 
nucleation,— hydrogen bonding,— insulator-conductor, sol-gel and glass transition a 3 i 4 i 5 i 6 i 7 i 8 i 9 
as well as bridging in granular materials^ are currently studied from this point of view. 
In all these cases, the system under study can be thought of as a collection of individuals 
(atoms, molecules, grains, etc.) that, with generality, we call particles. Most of the efforts 
have been based on lattice representations of the systems. The relative simplicity of lat- 
tice models allows for a wide variety of treatments, which extend from almost heuristic^ 
to quite rigorous.— Whatever the treatment is, the concept of connectivity between the 
particles plays an important role. 

Sometimes, however, a continuous description — where particles can occupy any point in 
a continuum phase space — is needed to reach a more realistic picture of the phenomena 
under consideration. For this context, the concept of connectivity has been generalized 
and adapted to describe clustering and percolation in continuum systems. The main ideas 
have been established in the pioneering works of HilU^ and Coniglio et a/.— Hill considers 
a partition of the whole system into subsystems of particles (the clusters) that satisfy some 
linking properties. The concept of cluster is thus directly related to the idea of bonded pairs. 
A bonded pair is a set of two particles that are linked by some direct mechanism. A cluster is 
then defined as a set of particles such that any pair of particles in the set is connected through 
a path of bonded pairs. We call these clusters "chemical clusters" to distinguish them from 
the non-pair-bonded clusters we have introduced in a previous worki^ — note, however, that 
this does not mean that clusters are necessarily formed through a true chemical bonding. 
A system is said to be in a percolated configuration if it contains a cluster that spans the 
system volume. 

From Hill's theory, we see that a connectivity criterion is needed in order to decide 
whether two particles are bonded or not. This connectivity criterion has to be defined in 
accordance with the phenomenon under study— i 16 ' 17 In the search for stable atomic clusters, 
which mark the onset of a phase transition in a monatomic gas, Hill proposed a simple 
energetic criterion (HE): two particles are bonded if their relative kinetic energy is less than 
the negative of their relative potential energy— In principle, this criterion takes into account 
the relative positions and velocities of the relevant pair of particles. For molecular fluids, 
instead of just their distance, the potential energy could in general depends on the direction 
and magnitude of the vector position of each of the two involved particles and their relative 
orientations. 

From a theoretical point of view, a criterion that involves the velocity of the particles 
prevents the straightforward integration of the momenta in the partition function, which is 
the great advantage of classical statistical mechanics. To avoid this obstacle, Hilli^ himself 
has proposed a velocity-averaged (VA) version of his criterion giving effective potentials 
between bound and unbound particles. 

The VA and the complete HE criteria have been used in computer simulations as well 
as in integral equations studies. For the VA criterion only the particle positions come into 
account, so it is suitable to both, Monte Carlo (MC) and molecular dynamics (MD) cal- 
culations. With respect to the integral equations approach, Coniglio et alM have obtained 
a connectedness Ornstein-Zernike (OZ) relationship for the pair connectedness function 
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gt(ri,r2) (see also a review from Stell).— This function is proportional to the joint prob- 
ability density of finding two particles belonging to the same cluster and at positions ri 
and r2, respectively. Therefore, by integrating 5^(1*1, r2), the mean cluster size S and the 
percolation density p p — i.e. the value of p for which S(p) diverges — can be obtained. Since 
Coniglio's theory deals only with the positions of the particles, the HE criterion cannot be 
implemented. Instead, the VA criterion was used by Coniglio et al^ to analytically calculate 
the percolation loci, for a potential made up of a hard core plus an attractive interaction, 
in a crude mean-field approximation. 

It is worth mentioning that most of the theoretical studie d 19 ' 20 ' 21 ' 22 on connectivity and 
percolation in continuum systems based on Coniglio's type equations were focused in the 
rather simple Stillinger's connectivity criterion.— This criterion states that two particles are 
bonded if they are separated by a distance shorter than a given connectivity distance d. In 
this case, d is an ad hoc parameter, which must be chosen on physical grounds. Although 
this criterion might be sensible in the study of certain insulator-conductor transitions, it is 
unrealistic regarding clustering in saturated vapors. 

A general theory which is appropriate for bonding criteria involving both, the momenta 
and positions of a pair of particles, has been developed by two of us.— The main object in 
our theory is the pair connectedness function g'(ri, r2, pi, P2) which is proportional to the 
joint probability density of finding two particles at positions ri and T2 with momenta pi and 
P2, respectively, and belonging to the same cluster. This function verifies also an OZ like 
relationship. In a previous paper— (thereafter denoted as I) we applied our general theory 
to study the complete HE criterion for the same model fluid considered by Coniglio et a/.— 
under the VA criterion. We also used the same simple closure relation proposed by Coniglio 
et al. More recently, we have reported 2 ^ the solution of our generalized connectedness OZ 
type relation for a Lennard- Jones fluid closed with a connectedness Percus-Yevick (PY) 
condition. We implemented a connectivity criterion which generalizes Stillinger's criterion 
in that a life time r for the bonds is required. 15 

In Ref. I we have also performed MD simulations of the Lennard- Jones fluid and have 
used both criteria (HE and VA) to define the clusters. We concluded that the VA criterion 
strongly overestimates percolation densities. We will partially revise these results here and 
will discuss some subtleties related to the identification of percolating clusters. Notice that 
MD simulations are convenient when the HE criterion is used to identify clusters since MC 
algorithms do not provide per se the velocities of the particles.— It should be mentioned 
that the HE criterion has been considered in MD studies of small clusters and the critical 
percolation behavior of Lennard- Jones fluids by several authors . 28 ' 29 ' 30 It has been suggested 
that the percolation line — the line that separates the temperature-density phase diagram 
into percolated and non-percolated states — might be experimentally observable . 30 ' 31 More- 
over, cluster analysis based on this criterion seems to be useful in locating the gas-liquid 
coexistence curve.— 

The main purpose of this work is to apply the generalized connectedness OZ type re- 
lationship closed with a connectedness PY condition for the Lennard- Jones fluid in order 
to obtain the pair connectedness function g\ iE (r 1 , r 2 , pi, P2) for HE clusters and thus the 
related cluster pair correlation function: 

9H E { r h r 2) = J P(ri, Pi)p(r 2 , p 2 )#H E (ri, r 2 , pi, p 2 )dpidp 2 , (1) 
where p(r 1 ,p 1 ) is N times the probability density of finding a particle at the phase space 
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configuration (ri, pi). 

We compare <7jj E (ri,r 2 ) so obtained with the function gl /A (r 1 ,r 2 ) for the VA criterion 
calculated using the integral equation that results when the OZ type relationship of Coniglio 
et aL lA is closed with a PY like conditions. Both functions — <?h E (i*i, r 2) and 5fy A (r 1; r 2 ) — are 
compared with the corresponding curves given by MD simulations. 

The paper is organized as follow. In Sec. II we present the model system and the two 
connectivity criteria, i.e. HE and VA, we will work on. Also we use this section to discuss 
some aspects about the MD simulations. The continuum clustering theories suitable for 
each criteria will be sketched in Sec. III. There, we briefly describe the integral equation 
for (yfy E (ri, r 2 , pi, P2) and its solution following Lado's orthogonal polynomials method . 26 ! 32 
Finally in Sec. IV we compare the theoretical results with those obtained from simulations. 
We then summarize and give the conclusions. 

II. MODEL SYSTEM AND ENERGETIC CRITERIA 

We consider a system of N particles whose configurations are given by their positions and 
momenta (rj,pj) (i = 1,...,N). The canonical (NVT) ensemble will be used throughout. 
We assume that particles interact via the Lennard- Jones pair potential 



where = |r^-| with = — iy 

The clustering criteria are expressed in terms of the bond conditional probability density 
P(r it j, Pij), say the probability density that two particles i and j are bonded under the 
condition that their positions and momenta are (rj,pj) and (rj,pj), respectively. 

A. HE criterion 

The original Hill's criterion (HE) identifies clusters by defining:— 



with Pi j the relative momentum: p it j = Pj — Pi- A maximum connectivity distance d has 
been added to the criterion in order to avoid unrealistic bonding. 

B. VA criterion 

By integrating the relative momenta weighted by the Maxwell distribution in the region 
where the relative kinetic energy is lesser than minus the pair potential, the momenta are 
eliminated and the VA bond conditional probability density is obtained : 13 ^ 4 





Pre^zj, Pi,j) — 



1 pfj/4m < — t>(rj, ij) and r^- < d 
pfj/4m ^ — v(ri,Tj) or r it j > d 



(3) 




(4) 



where T[a] is the gamma function and 7 [a, x] is the incomplete gamma function. 
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C. Molecular dynamics calculations 



We consider a system of Lennard- Jones particles in a cubic box with periodic boundary 
conditions in the NVT ensemble and use a leap-frog algorithm with velocity correction.— 
The time step is chosen as At* = Ata^ 1 ^JkgT / (em) = 0.01. Quantities are averaged over 
10 3 configurations chosen every 100 At after equilibration. A cut off distance equal to 2.5a 
was used in the pair potential. 

In the VA case, for each pair of particles that satisfies < and ry < d, we generate 

a random number z, between and 1. If z < 7[3/2, — t>(rjj)//csT]/r[3/2] we consider that 
the particles form a bonded-pair, otherwise we do not. Note that this criterion can also be 
used in MC simulations because it does not require information on the particle velocities. 
To identify the clusters from the list of bonded-pairs we use the Stoddard's algorithm . 33,3 " 

A system is said to be in a percolated state if a cluster that spans the replicas is present 50 
percent of the time.— It is known that this criterion yields results that are marginally affected 
by finite size effects.— Then, a percolation transition curve, which separates the percolated 
from the non-percolated states of the system, can be drawn above the coexistence curve in 
the T — p phase diagram. 

In Fig. 1, the percolation loci for HE and VA connectivity criteria are presented for 
d = 3a. These simulations were performed with iV = 1372 particles. The gas-liquid 
coexistence curve obtained by Panagiotopoulos 37 and the MC liquid-solid coexistence curve 
of Hansen and Verlet 3 ^ are also shown. The MD results for the HE criterion are similar to 
those obtained by Campi et al. 30 These authors consider that the system is on the percolation 
line if the second moment of the cluster size distribution that excludes the largest cluster 
n'(s) reaches its maximum. In Fig. 2 we show the dependence of the calculated percolation 
density with system size and connectivity distance d. Extrapolation to infinite systems can 
be obtained by fitting a straight-line to a plot of p p versus L~ 1//l/ .— We have used the universal 
value of v = 0.88 ± 0.02 reported by Gaunt and Sykes 39 for three-dimensional systems. The 
error due to finite size effects in the calculated value of p p for the 1372-particle system is of 
1.0%. This correction is smaller than the size of our symbols in Fig. 1. Also from Fig. 2, 
we see that the effect of the connectivity distance d on p p is negligible for d > 2.5. 

As we can see from Fig. 1, the VA criterion is a very good approximation to the full HE 
criterion as far as the percolation loci is concerned. In Ref. I, we reported a VA percolation 
line that was located at much larger densities, and concluded that the approximation was 
rather poor. The revised results reported here show that this is not the case. The source 
of the error in Ref. I comes from the way percolating clusters are detected according to the 
Seaton-Glandt prescription. 35 All clusters in a given configuration are first identified by the 
Stoddard's algorithm, then each cluster is analyzed separately to determine if its replicas 
are connected with one another. Since the VA criterion implies the use of random numbers 
to decide whether two particles are connected, the second step where each separated cluster 
is analyzed for percolation needs to reuse the same random numbers generated when it was 
first identified. This subtlety was overseen in Ref. I, which led to the incorrect identification 
of actual percolating clusters as disconnected replicas. 

III. CLUSTER PAIR CORRELATION FUNCTIONS 

In the remainder of the paper we restrict our attention to the cluster pair correlations for 
the HE and VA energetic criteria. We calculate them by using the above mentioned integral 
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equations and MD simulations. Thus, this section will be devoted to pose the integral 
equations for the cluster correlation functions <7H E (ri, r 2 ) and g VA (ri, r 2 ) and briefly discuss 
their solutions. 

A. VA criterion 

In order to study clustering in a system composed of N classical particles interacting via 
a pair potential t>(ri,r 2 ), Hill separated the Boltzmann factor e(ri,r 2 ) = exp[— /3v(ri, r 2 )], 
into bonded (t) and unbounded (*) terms: 13 e(r 1; r 2 ) = e^(r!,r 2 ) + e*(ri,r 2 ). As usual (3 = 
1/kgT, being kg the Boltzmann constant. Since e^(r 1; r 2 ) represents the basic probability 
density of finding two particles bonded and at positions r 1 and r 2 , this separation yields 
a diagrammatic expansion for the partition function in terms of "chemical" clusters. We 
express Hill's separation as follows 

e\r 1 , r 2 ) = F(r x , r 2 ) exp[-f3v(r h r 2 )] (5) 

e*(r 1; r 2 ) = [1 - P(r 1; r 2 )] exp[-/3i;(n, r 2 )] (6) 

where P(r 1; r 2 ) = -PvA( r i,2) is given by Eq. (jlj) in the case of the VA energetic criterion. 

Fugacity and density expansions have been found, within Hill's formalism, by Coniglio 
and co-workers^ for the pair connectedness function ^(ri, r 2 ) = <7 VA (ri,r 2 ). As it was 
already mentioned, this function is proportional to the joint probability density of finding 
two particles belonging to the same cluster and at positions rj and r 2 , respectively. Moreover, 
by collecting nodal and non-nodal diagrams in these expansions an OZ-type relationship is 
obtained 

g j (r 1 ,v 2 ) = c f (ri,r 2 ) + p J c j (v 1 ,r 3 )g f (v 3 ,r 2 )dr 3 . (7) 

The function c^(r 1; r 2 ) is the direct pair connectedness function. By posing a closure 
relation, an integral equation for g*(ri,r 2 ) is obtained. Here, we use the more reliable 
closure available, i.e. the PY-like relation 5 

g\r 1>2 ) = [f*(r 1)2 ) + l]^(ri, 2 ) - c+(r 1)2 )] + exp[^(r li2 )]^(r 1|2 )/t( rii2 ). ( 8 ) 

In Eq. (jHJ), /*(ri )2 ) = e*(ri j2 ) — 1 = exp[— j3v (ri j2 )][l — -Pva(^i,2)] — 1 is the unbound Mayer 
function and g{r\ )2 ) is the thermal pair distribution function (PDF). 

In order to solve the integral equation given by Eqs. (|7j) and (jHJ), for the Lennard- Jones 
potential we have implemented Labik's numerical algorithm.— 

B. HE criterion 

1. The integral equation 

We summarize here the basic theory that we have developed^ to describe the clustering 
and percolation for clusters whose bond definition depends on the positions and momenta 
of the two particles under consideration. 
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For a system of N classical particles that interact through a pair potential v(vi,rj), we 
define a density correlation function p(r\, r 2 , pi, P2) that is N(N — 1) times the probability 
density of finding two particles at the phase space configurations (ri, pi) and (r 2 , p 2 ) 
respectively: 



P(ri,r 2 ,pi,p 2 ) " 



h 3N N\Q(N,V,T) 

N 2 



X 



/ n ex p[-^^] n n «phm*« r^d^-v- 2 . (9) 



i=l i=l j>i 



Here h is Planck's constant and Q(N, V,T) the canonical partition function of the system. 
Then, in the same spirit of Hill and Coniglio et al., 5,13 we separate exp[— (3v (r^, r^)] into 
connecting and blocking parts, 

exp[-(3v (r h rj)] = f(r h r j; p { , pj) + f*(r h r j; p l; pj) + 1. (10) 

Here /'(rj, Vj, Pi, pj) represents the basic probability density that two particles in con- 
figuration (rj, Tj, pi, pj) are bonded. We will sometimes use the shorthand notation 
/ 7 (rj, rj, pi, pj) = f?j, where 7 can be either f or *. Substitution of Eq. (fTUj) in Eq. 
© yields 

N(N-l) 

p(ri ' r2 ' Pl ' P2) = h^N\Q(N,V,T) eXp[ -^ r2)] 



c N 2 

x / n ^-^fy E^n fhf:^ N - 2 ^- 2 ' (ii) 
i=i ~ in 



where the sum is carried out over all possible arrangements of products of functions f\j and 

We note that the functions • and f*j can depend on the momenta as well as on the 

positions of the two particles, but the sum of fjj and f*j must be momentum independent 

in order to conform to Eq. (J1(J)) . Except for this last condition, the functions /|- and f*j 
are otherwise arbitrary for thermodynamic purposes. Of course, we choose them in such a 
way that the desired definition of bonded particles for HE clusters is achieved, i.e., 

flj = e W [-i3v(r i; j)}P(v i) j,p i;j ) (12) 



f*j = exp[-/Mr„)][l - P(n,j, Pij)} ~ 1 (13) 

where P(x ij j,p ii j) = P}m(*i,j, Pi,j) is given in Eq. fl3J). 

Each term in the integrand of Eq. (fTTj) can be represented as a diagram consisting of 
two white ei and e 2 points, N — 2 black points and some • and f*j connections except 

between the white points. Here we take e; = exp[— /Jg-M. White points are not integrated 
over whereas black points are integrated over both their positions and momenta. All the 
machinery normally used to handle standard diagrams in classical liquid theory^- can now 
be extended to treat these new type of diagrams. By following Coniglio's recipe to separate 
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connecting and blocking parts in the PDF, g(ri,r 2 ) = g^(ri, r 2 , pi, p 2 ) + g*(r\, r 2 , Pi, P2), 
we obtain an OZ-like integral equation for g^(ri, r 2 , pi, p 2 ), 

<7 f (ri, r 2 , pi, p 2 ) = c f (ri, r 2 , p a , p 2 ) 

+ y p(r 3 ,P3)c t (ri,r3,pi,p 3 )5( t (r3,r 2 ,p3,p 2 )rfr3rfp3. (14) 

Here p(ri, pi)p(r 2 , p 2 )g^(ri, r 2 , pi, p 2 ) is N(N — 1) times the joint probability density of 
finding two particles at positions ri and r 2 with momenta pi and p 2 , respectively, and 
belonging to the same cluster, where the bonding criterion is given by Eqs. (jl2j) . ()13|) and 
©, while 

p(ri,Pi) = jfZTJ J P( r i> r 2,Pi,P2)dr 2 dp 2 . (15) 

The function c^(r 1; r 2 , p 1; p 2 ) denotes the sum of all the non-nodal diagrams in the dia- 
grammatic expansion of g'(rx, r 2 , Pi, p 2 ). We recall here that a nodal diagram contains at 
least one black point through which all paths between the two white points pass. For a 
homogeneous system, we have 



£ f (ri2, Pi, P2) = c f (r 12 , Pl , p 2 ) + 



{2^mk B Tf/ 2 

f v 2 

x / exp[-/3^]c t (r 13 ,p 1 ,p3)^ t (r3 2 ,p3,p 2 )c/r 3 dp3. (16) 

To obtain a closed integral equation with Eq. (fTlj) or Eq. (flljj) . we need a closure 
relation between gt(ri, r 2 , pi, p 2 ) and c^(ri, r 2 , pi, p 2 ). Here we will use the PY approxi- 
mation g(ri, r 2 ) exp[/?v(ri, r 2 )] = 1 + iV(ri,r 2 ), where the function iV(ri,r 2 ) is the sum of 
the nodal diagrams in the expansion of g(ri,r 2 ). Separation into connecting and blocking 
parts, ^(r!, r 2 ) = gi(n, r 2 , p 1; p 2 ) + p*(n, r 2 , Pi, p 2 ) and jV(ri,r 2 ) = iV^ri, r 2 , p 1; p 2 ) + 
iV*(ri,r 2 ,p a ,p 2 ), yields 

flf^ri, r 2 , p 1; p 2 ) = [f*(r u r 2 , p ls p 2 ) + 1] [^(r^ r 2 , p 1; p 2 ) - c f (r 1; r 2 , p 1; p 2 )] 

+ exp[/3w(ri, r 2 )]g(r 1 , r 2 )/ t (ri, r 2 , p 1} p 2 ), (17) 

or, for a homogeneous system, 

# f (ri 2) Pi, Pa) = [/*(ria, Pi, P2) + 1] [# f (ri 2 , pi, p 2 ) - c f (ri 2 , pi, p 2 )] 

+ exp[^(ri 2 ]5((ri 2 )/ t (ri 2 , p h p 2 ). (18) 

Equation ([14]) joined with Eq. (fT7j) or Eq. (fTB|) joined with Eq. (JTHJ) give a closed set of 
equations for g Jf (r 1 , r 2 , p x , p 2 ). 

From the function g^ E (ri, r 2 , pi, p 2 ) = g^(ri, r 2 , pi, p 2 ) we define the pair correlation 
function for energetic clusters <7h E (i"i, r 2 ) according to Eq. (JTJ). 
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2. Solution of the integral equation 



Equivalence with an integral equation for polarizable fluids Our problem con- 
sists in solving Eq. (j!6|) for <^(r 12 , pi, P2) closed by the connectedness PY relation (fTS|) with 
/^(rj, Tj, Pi, pj) and /*(rj, r^, pj, pj) given by Eqs. (|T2|) and (fT3|) . In the closure relation 
(|18|) . (?(r 12 ) is the thermal PDF of the system. In this work we take g(r 12 ) from the solution 
of the thermal OZ equation in the PY approximation— 

An equation mathematically equivalent to Eq. ()16j) has been previously solved by Lade 32 
in the study of nonpolar polarizable molecules. Explicitly, the equation considered there, 
which is a generalized OZ equation, relates the fluid total correlation function (TCF) 
h(ri2, pi, p 2 ) = <?(ri2, pi, p 2 ) — 1 (with g(ri 2 , pi, P2) the PDF) and the direct correlation 
function (DCF) c(ri 2 , pi, p 2 ), 



where Pi denotes the instantaneous dipolar moment induced on molecule i by the remaining 
molecules of the system. The function / (p) gives the instantaneous dipolar moment thermal 
distribution which is assumed to have a Gaussian form 



where a is the effective polarizability of the molecules. 

We observe that Eqs. (fTB|) and (|T9"j) are the same if we identify h with g\ c with c', the 
induced dipolar moment p; with the kinetic momentum pi and the polarizability a with the 
particle mass m. There are, however, some differences between the connectivity problem 
and the polarizable- molecule problem. The form of f(p) does not need to be Gaussian for 
polarizable molecules; moreover, f(p) is coupled to the TCF. Therefore, the value of the 
effective polarizability a depends on the density and temperature of the system. In the 
connectivity problem, however, the equivalent of f(p), p(r,p)/p, is intrinsically Gaussian 
and independent of the thermodynamic macrostate of the system. 

Another difference between the connectivity problem here and the problem described by 
Lado is that our closure relation must be complemented with the condition given by Eqs. 
()12j) and (|13|). In addition, the closures are different. Here we consider the connectedness 
version of PY whereas an almost exact relation between DCF and TCF (van Leeuwen- 
Groeneveld-De Boer— exact relation with approximate bridge function) is used by Lado.— 
Nevertheless, these differences do not affect the general method of solution developed by 
Lado and we can apply the same principle of expansions in orthogonal functions. 

Thus, following Lado , 26 ' 32 we start by reassigning the unknown function to be the indirect 
correlation function 



M r 12,Pl,P2) 



c(ri2,Pi,p 2 ) 




(19) 




7 (ri2, Pi, P2) = g\r n , Pi, P2) - c j (ri2, pi, p 2 ), 
rather than g^(r 12 , p 1; p 2 ), and rewriting Eq. ()16|) in Fourier representation, 



(20) 
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f (k ' Pl ' P2) = (2nmk B T)W I 

[f(k, Pl , p 3 ) + ^(k, Pl , p 3 )] c t (k, p 3 , p 2 ). (21) 

The closure given by the PY relation [Eq. (|18|) ] together with the conditions (|12j) . (|13|) 
and Q yield 



c t (ri2,Pi,p 2 ) 



g(r 12 ) - 7 t (ri2,Pi,P2) P?, 2 / 4m < ~v(r 12 ) and r 12 < d 

(exp{-(3v(r 12 )] - 1) 7 t (r 12 , p x , p 2 ) P^/4m > -v(r 12 ) or r lj2 > d 

(22) 

The connectivity part of the PDF is then computed from y as 



, „ n] J 9(ru) Pi, 2 / 4m < ~v(r 12 ) and r 12 < d 

9 lr 12 , Pl ,p 2 j - <; exp[ _^ (ri2)]7 t( ri2)Pl)P2 ) p 2. /4m > _ u(ria) or ri2 > d (-0) 

The Fourier transform in Eq. (|2*T|) and its inverse are defined as 

/»= [ drf(r)e-^, (24) 



/ (r) = —-3 J dkf (k) e* kr . (25) 

The standard method for solving Eqs. (j2~Tj) and (jUj) is to explicitly break out the angular 
dependence of all functions in the form of expansions in spherical harmonics.— 

Expansion of the pair functions in orthogonal polynomials The essential point 
in the integral equation solution method^ is the expansion of all the pair functions, like 
7*( r i2> Pi; P2), in terms of orthogonal polynomials. First we expand 



7 f (ri2,Pi,P2) = 7 f (r,pi,p 2 ,u;i,u; 2 ) 

= 4vr l\ lh m,{ r ,Vi,V2)Y hm {ui)Y hm (u 2 ) , (26) 

where uj\ and u 2 are the directions of the momenta pi and p 2 , m = — m, and m = —I, —I + 
1, In this and similar expressions, the vector r 12 has been implicitly chosen as the z 
direction in the specification of the Euler angles u = (8, <fi). The spherical harmonics satisfy 
the orthogonality condition 

j duY lm tu)Yf rj ,{u)=8 ll ,5 mm ,i (27) 
so that the coefficients of the expansion (J26|) are immediately obtainable as 

lhhm{ r iP^P2) = 4^ / dhJ 1 du 2 ^\r,pi,p 2 ,u u u 2 )Y hm {u x )Y^ m {u 2 ) . (28) 
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Similarly, we can break out the kinetic momentum dependence in the form of expansions 
in polynomials of p, 



lihrni^PliPz) = E Ihhm ( r ) Qnih (Pi) Qn 2 l 2 (p 2 ) 



ni,ra 2 



which are constructed to be orthogonal with Gaussian weight function 

1 



f(p) 



(2vrm//3) 3 / 2 



exp[— (3p /2m] 



(29) 



(30) 



namely, 



00 

4vr J dpp 2 f (p) Q nl (p) g n ;, (p) = 5 r 



(31) 



The coefficients of the expansion are then again obtainable by quadratures, 



t nin.2 
1 (1/2T1 



rfpirfp 2 [^p\f (pi)] [47rp 2 / (p 2 
x 7Ut»( r >Pi>P2)GmJi (pi) Q„ 2 ; 2 (p 2 ) • 



(32) 



Given the Gaussian form of the weight function / (p), the associated polynomials are^i 

r(i(n-0 + i)r(|) 



Qnl (P) 



r(|(n + + f) 



1/2 'fl^x'/a 

PP ^ ^2+1/2 / PP 



2m 



J L (n-/)/2 



2m 



(33) 



where (£) are the associated Laguerre polynomials^ and T (z) is the gamma function. 
Accordingly, all the functions in r-space are expanded in the form 

7 f (r, pi, pa) = 4tt 7 J i ™ 1 ™ 2 (r) Q nih (p x ) Q n2h (p 2 ) Y, im (wi) Y; 2 ™ M , (34) 

ni, n 2 . ii, i2i wi 

where the 2 axis is along r and the summation indices satisfy the constraints 



n = 0,l,2,..., 



I = n, n — 2, n — 4, 1 or 0, 
m = 0, ±1, ±2, ±Z. 



The coefficients of Eq. (13411 can be obtained as 



(35) 



lllZ 2 (r) = 4tt J d Pl dp 2 f (pO / (p 2 ) 7 t(r, Pl , p 2 ) 
X Qnih (Pi) Qmh (Pa) (^l) ^'m M 



(36) 
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with / (p) given by Eq. (|30|). The complete orthonormality condition is 



4vr j dp/ (p) (p) g nY ( P ) Y lm ( W ) y; m , (w) = <UMw • ( 37 ) 

The functions in k can be expanded in a similar way. Setting the z axis along k, we write 



7 f (k, pi, pa) = 4vr TijSJT (Pi) Qnafa (P2) M ^ 2 m (w 2 ) 



(38) 



ni, n 2 , ii, ia, m 



However, the angles u% and co> 2 are referred to different axes in Eqs. (|3~4*j) and (jHHJ), so that 
the coefficients in these expansions are not themselves mutual Fourier transforms. 

Introducing the expansion for 7'(k, pi,p 2 ) and the corresponding expansion for 
ct(k, pi, p 2 ), one finds that the OZ-like equation in Fourier space [Eq. (J2IJ)] goes over 
into a set of matrix equations for the respective coefficients, 



7&r w 



n3,h 



~t ™3™2 
C l 3 hm 



(k). 



(39) 



Numerical procedure To obtain a numerical solution for the set of equations (jlfij) 
and (jTSjl one needs the discrete versions of the expansion for 7*(r, Pi, P2) [Eq. ()34|l] and the 



quadratures for the coefficients 7] 



t n l n 2 



[Eq. (|3*Sj)]: these are 



and 



7 t (r,z 1 ,i 2 , fci,fc 2 , j) = 4vr 7 / t i ™ 1 ^ 2 (r) Q niJl Q„ 2 j 2 (i 2 ) 

ni, rt2, ii, ^2, m 

x P, im (fei) P, aW (fc 2 ) v m T m (7) (40) 



~t una /y 



2 iw (n) u> (z 2 ) u> (fei) to (fc 2 ) iw (j) 7 f (r, ii, z 2 , fci, & 2 , j) 

il,i 2 , fci,fca,i=l 

x Q mJl (ii) Q na , 2 (i 2 ) 7\™ (A*) TV (A&) (-l) m T m (j) . (41) 

In Eq. f)4(J|l . u = 1 and i/ m = 2 for m > 0. In Eq. (jHJ, Gaussian quadratures are being used, 

with the argument i standing for U = j3pf/2m, the zth root of L]( 2 (t), k for Xk = cos 9k, the 

1/2 

kth root of Pn v (x), and j for yj = cos <pj, the jth root of T^ p (y), where (t), Pn v (x) , and 
T/v (y) are the associated Laguerre, Legendre, and Chebyshev polynomials, respectively, all 
of order N p ; here the associated Legendre functions Vi m {x) are normalized to 2. The w are 
the corresponding Gaussian weights, 



w (1) 



U 



rl/2' 



-1 



W 1 



(*)= (i-x 2 fe ) p;> fc ) 



(42) 
(43) 
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w(j) = N-\ (44) 

where the prime denotes derivative. 

The solution follows an iterative procedure. The preparatory stages of the calculation 
consist of (i) computing the thermal PDF g (r 12 ) for the Lennard- Jones fluid over a suitable 
mesh using the PY equation, (ii) reducing the momentum space to the discrete set of points 
Pi,k,j = (pi, 9k, 4>j) with i,k,j = 1,2, ...,N p , and (iii) identifying the subset of states - 
within all possible configurational states (r 12 , pi, p 2 ) of a pair of particles — that correspond 
to a bonded pair. 

We construct a logical array B(ri 2 , Pi ; i,fcj, P2;i,fcj) of dimension seven whose value is 
TRUE if the configurational state of the pair of particles corresponds to a bonded state, 
i.e., if pi 2 /4m < — u(ri 2 ) and ri 2 < d. If instead this condition is not satisfied, then 
B(r 12 ,Pi;i,fc J j,P2;i ) fe,j) is FALSE. 

The iterative solution of Eqs. ()39|) and (|22j) starts by guessing the initial values of the 
coefficients 7y^ 2 (^12)- Then, if B(r 12 , Pi ; i,fcj, P2;i,fcj) is TRUE, following Eq. we take 

tmna/ \ / ^( r i2) if n 1 = n 2 = l 1 =l 2 = m = 
^™ (ri2) = \0 otherwise. (45) 

If instead B(ri 2 , Pi ; i,jfe,j, P2;i,fe,i) is FALSE then, following Eq. (|2*3*|l. we take 

ff& fa) = exp[-^(r 12 )] 7 ;;- 2 (r 12 ) . (46) 

Knowing 0/ (^12) and ^i\^m ( r i2) f° r an the mesh points and allowed indices, we can 
calculate [see Eqs. (|20|) or (122J) ] 

<r 2 r (-12) = 9{IZ 2 (ru) ~ lUT (r u ) • (47) 

We now need to transform the coefficients c]^ 1 ^ 2 (r 12 ) in real space into coefficients 

^ Tm 2 (^) * n F° ur i er space. However, as we have mentioned, they are not themselves Fourier 
transforms of each other. Thus, we have to assemble the complete function first using the 
equation analogous to (j4*Uj) for c^(r, i±, i 2 , k±, k 2 ,j) and then use a generalized fast-transform 
algorithm 32 to calculate $(k, i±, i 2 , ki, k 2 ,j). Using the equation analogous to (|4*Tj) in k-space 
we then have the coefficients c^m 2 (^) f° r the complete set of indices and all the values of 
k on an adequate mesh. The coefficients 7J f 1 ™ 2 (k) are then easily calculated by using the 
OZ-like equation in Fourier space [see Eq. ()39|)]. Again, we assemble the complete function 
j* (k, ii, i 2 , ki, fc 2 , j) [using the Fourier space version of Eq. (|40j)]. The inverse transform 
7H r i2) h, h, ki, k 2 ,j) is calculated with the fast-transform algorithm and so new coefficients 
TiJ*m ( ri2 ) [obtained from Eq. (|4*T|)] are available to start again the iterative cycle. The 
iterations end when convergence is reached, as measured by 



J (s+l)tn iteration L 



t "1™2 



- sth iteration 



< e (48) 



for the complete set of indices. The tolerance e is set to 0.0001 . 

The pair correlation function for an energetic cluster [see Eq. (JIJ] is finally given by 

#H E ( r i2) = 9ooo ( r i2) , (49) 
where the orthonormality condition [see Eq. 1)370] has been used. 
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IV. RESULTS 



Firstly, as a complement to Fig. 1, we show in Fig. 3 the cluster pair correlation functions 
5^(7-12) and # VA (r 12 ) obtained from MD using N = 4000 particles for (p*,T*) = (0.24,1.4), 
(p*,T*) = (0.42,1.4) and (p*,T*) = (0.429,1.4) where p* = pa 3 and T* = k B T/e. The two 
last points correspond, respectively, to the VA and HE percolation loci for T* = 1.4. The 
main peak in the cluster correlation functions is higher for the VA criterion than for the HE 
criterion. This implies that there is a larger tendency to consider as directly connected two 
neighbor particles by this criterion. However, as it can be appreciated in Fig. 3, gy A (ri 2 ) 
falls faster than g^^r^) for larger r. This contrasting behavior can be understood by 
analyzing the cluster size distribution function n(s), which gives the number of clusters in 
the system consisting of s particles. Figure 4 shows n(s) for the same three state points as 
Fig. 3. We can see here that the VA criterion identifies a larger amount of clusters than the 
HE criterion up to a certain size — which depends on the density. However the HE criterion 
always identifies some clusters larger than the largest clusters identified by the VA criterion. 
This means that the cluster correlation function is more long ranged for the HE clusters. 

Figs. 5 and 6 show the theoretical cluster correlation functions g\^{r 12) and gv A (ri 2 ), 
calculated by solving the corresponding integral equations [Eqs. ()16p .(fT% |) and (0),©, re- 
spectively] following the methods indicated in the previous section, together with the corre- 
sponding MD simulation results. We show the curves at temperatures T* = 1.4 and 3.0 and 
densities p* = 0.24 and 0.55, respectively. This values are rather far from the percolation 
loci since we have been unable so far to obtain convergence of the numerical algorithms at 
higher densities when the HE criterion is considered. 

With generality, Figs. 5 and 6 say that the theoretical results follow quite well the trends 
of the corresponding simulations for each criterion. 

V. CONCLUSIONS 

We have shown that, contrary to our previous results (Ref. I), the VA energetic criterion 
is, in general, a good approximation to the full HE criterion for estimating the percola- 
tion loci in a Lennard- Jones fluid. However, the cluster correlation functions are somewhat 
different in the VA case. We have obtained the cluster pair correlation functions for both 
energetic criteria through the numerical integration of connectedness OZ integral equations. 
In particular, we have used a generalization of the integral equations that allows the imple- 
mentation of the HE criterion. The theoretical results agree rather well with the simulations. 
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Figure Captions 



Figure 1. Coexistence and percolation curves for the Lennard- Jones fluid. Open dia- 
monds correspond to the gas-liquid coexistence.— Open squares correspond to the fluid-solid 
coexistence.— The percolation loci for the HE criterion (open circles) and for the VA crite- 
rion (open triangles) are compared. Solid circles correspond to the percolation loci for the 
HE criterion from Campi et alM Lines are only to guide the eye. 

Figure 2. Percolation density for the HE criterion as a function of L~ x l v . L is the 
simulation box length and v = 0.88.— The largest system corresponds to 4000 particles. 
The system size used for the calculation of the percolation curves in Fig. 1 is indicated by 
an arrow. The inset shows the percolation threshold as a function of d for the 1372-particle 
system. The arrow shows the value of d used in the rest of the paper. 

Figure 3. Connectedness correlation function for the HE criterion (solid line) and for 
the VA criterion (dotted line) at T* = 1.4 for various densities. The results for p* = 0.42 
and 0.429 correspond to the percolation loci for the VA and the HE criterion, respectively. 

Figure 4. Cluster size distribution for the HE criterion (black solid line) and for the VA 
criterion (red dotted line) at T* = 1.4 for the same densities as in Fig. 2. 

Figure 5. Connectedness correlation function for the HE criterion (solid line and open 
circles) and for the VA criterion (dotted line and open triangles) at T* = 1.4 and p* = 
0.24. Lines correspond to the solution of the connectedness OZ equations closed with the 
connectedness PY relation. Symbols correspond to MD. 

Figure 6. Connectedness correlation function for the HE criterion (solid line and open 
circles) and for the VA criterion (dotted line and open triangles) at T* = 3.0 and p* = 0.55. 
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